El modelado de distribución de especies es una herramienta clave en ecología y biogeografía que permite predecir y entender cómo y por qué ciertas especies se distribuyen en un hábitat o espacio geográfico determinado. Tradicionalmente, utiliza datos demográficos de las especies (bien presencia/ausencia, bien abundancia) junto con variables ambientales para construir modelos que simulan su distribución en diferentes escenarios. Estos modelos son fundamentales para la conservación de la biodiversidad, el manejo de ecosistemas, la restauración de hábitats y el estudio de los impactos del cambio climático sobre las especies. Además, proporcionan información crucial para la toma de decisiones en, por ejemplo, la planificación ambiental y el diseño de áreas protegidas.
El problema de la cuarta esquina (fourth-corner) trata de relacionar los rasgos de los organismos con el entorno en el que estos se encuentran. Es una aproximación más holística al modelado de distribución de especies, que clásicamente se realiza teniendo en cuenta la relación directa de las especies con el ambiente. ¿Por qué recibe el nombre de la cuarta esquina? Se debe a que involucra tres matrices para calcular una cuarta, como vemos en el siguiente esquema:
Y, ¿por qué decimos que es un problema? Esto se debe a que no es nada sencillo calcular dicha matriz. En los últimos 30 años ha habido dos propuestas principales:
Aunque ambas soluciones cumplen diferentes propósitos, no proporcionan una solución cuantitativa para la cuarta esquina y, por tanto, no nos permiten entender la magnitud de las interacciones entre rasgos y ambiente. En los pasos siguientes, veremos cómo podemos hacer esto utilizando modelos lineales generalizados con las tres matrices, siguiendo la solución propuesta por Brown et al. (2014).
En primer lugar, cargaremos las librerías de R que vamos a utilizar a
lo largo de nuestros análisis. Las librerías o paquetes
son recopilaciones de funciones, generalmente sobre un temática
específica, que nos sirven para realizar procesos en R, como manejo o
análisis de datos, de una forma estandarizada y rápida. Aunque estos
paquetes suelen estar escritos por personas que se dedican a la
programación o a la ciencia, R es una comunidad de software libre, así
que cualquier persona que tenga una metodología y quiera compartirla con
la comunidad de usuarios de R puede crear su propio paquete. Los
paquetes de R pueden instalarse directamente a través del CRAN utilizando la función
install_packages() con el nombre del paquete como argumento
de la función, entre comillas (por ejemplo:
install.packages("tidyverse")). En ocasiones, también
pueden instalarse directamente desde las cuentas de GitHub de los usuarios. Los paquetes solo
necesitan instalarse una vez en un ordenador (aunque conviene estar
informado sobre posibles actualizaciones). Sin embargo, cada vez que
iniciamos una nueva sesión de R hemos de cargar esos paquetes. Esto se
hace usando la función library(), que es la forma de
decirle a nuestro programa que queremos utilizar las funciones que el
paquete contiene.
A continuación, pues, instalamos las librerías a utilizar:
install.packages("tidyverse") # manipulación, transformación y dataviz
install.packages("mvabund") # manupilación y análisis de datos multivariantesY las cargamos:
# fíjate en que ya no hacen falta las comillas para cargar los paquetes
library(tidyverse)
library(mvabund)Estos paquetes son:
tidyverse (Wickham et al. 2019), que es
en realidad una colección de paquetes de programación funcional que
sirven para la manipulación de datos de forma más eficiente, utilizando
tuberías (pipes) para enlazar elementos. La pipe es
%>%, que veremos mucho en el código en este documento.
Aunque parezca poco tostón escribir %, > y
%, podéis insertarla automáticamente presionando las teclas
shift + control + m en Windows (o
shift + command + m en Mac).
mvabund (Wang et al. 2012), que es un
paquete para modelizar datos de abundancia en ecología de comunidades y
que contiene las funciones que utilizaremos en el análisis de
fourth-corner.
Vamos a crear matrices con datos ficticios (lo que se llama en inglés datos dummy) para aprender cómo se utiliza esta técnica. Lo primero que haremos será seleccionar cuántas especies, cuántos rasgos funcionales, cuántas variables ambientales y cuántos sitios o localidades queremos:
species <- 20 # número de especies para nuestros datos dummy
sites <- 10 # número de sitios o localidades para nuestros datos dummy
traits <- 15 # número de rasgos funcionales para nuestros datos dummy
env_var <- 15 # número de variables ambientales para nuestros datos dummyY, a continuación, vamos a sembrar una semilla (seed). ¿Qué
quiere decir esto? La semilla es un número natural que transforma lo
estocástico en determinista. Esto es, al trabajar con funciones que
generan datos aleatorios (por ejemplo, runif(n), siendo
n el número de elementos a generar), si no hemos plantado
una semilla, cada vez que corramos la función los resultados generados
cambiarán (¡probad a correr runif(10) sin semilla varias
veces!). Sin embargo, cuando plantamos una semilla usando la función
set.seed(), el resultado generado será siempre el mismo, y
sólo cambiará si cambiamos la semilla. Aquí utilizaremos el número 5558
como nuestra semilla por defecto.
Finalmente, vamos a crear las matrices necesarias para los análisis. Tenemos en mente que, aunque en este ejercicio estamos generando datos ficticios, en la práctica estas matrices comprenderían datos reales de nuestro sistema de estudio. La cuestión es que estos datos deben estar organizados en las siguientes matrices con la siguiente estructura:
Matriz de abundancia de especies por localidad (L), con 20 especies en columnas y 10 localidades en filas.
# crear una base de datos de combinaciones de especie y sitio
L <- expand.grid(species = paste0("sp", 1:species),
sites = paste0("site", 1:sites))
# añadir una columna con los datos de abundancia
L$value <- rpois(nrow(L), 30)
# incorporar los datos a una matriz
L <- L %>%
spread(key = species, value = value) %>%
column_to_rownames(var = "sites") %>%
as.data.frame()
# observar las primeras filas de dicha matriz
head(L)Matriz con datos de rasgos funcionales por especie (Q), con 15 rasgos en filas y 20 especies en columnas.
# crear una base de datos de combinaciones de rasgos y especies
Q <- expand.grid(species = paste0("sp", 1:species),
traits = paste0("t", 1:traits))
# añadir una columna con los datos cuantitativos
Q$value <- round(runif(n = nrow(Q), 0, 1), 3)
# incorporar los datos a una matriz
Q <- Q %>%
spread(key = species, value = value) %>%
column_to_rownames(var = "traits") %>%
as.data.frame()
# observar las primeras filas de dicha matriz
head(Q)Matriz con datos de variables ambientales por localidad (R), con 15 variables ambientales y 10 localidades.
# crear una base de datos de comb. de variables ambientales y sitios
R <- expand.grid(environ = paste0("env", 1:env_var),
sites = paste0("site", 1:sites))
# Add a column with random quantitative data from 0 to 1
R$value <- round(runif(n = nrow(R), 0, 1), 3)
# incorporar los datos a una matriz
R <- R %>%
spread(key = environ, value = value) %>%
column_to_rownames(var = "sites") %>%
as.data.frame()
# observar las primeras filas de dicha matriz
head(R)El objetivo de esta sesión no es aprender R, ni mucho menos dominar la visualización de datos (que requiere mucha práctica), pero sí me gustaría que os lleváseis una idea de que estas herramientas están disponibles y no son tan innaccesibles como pueden parecer.
Lo primero que vamos a hacer es observar cómo es la distribución de
la abundancia de organismos en las diferentes poblaciones. Para ello,
vamos a utilizar el paquete ggplot2, perteneciente al
tidyverse, del que ya hablamos anteriormente. La función
ggplot() es una de las más poderosas del lenguaje R para
crear visualizaciones de datos. Funciona por capas: primero,
ggplot crea un espacio o lienzo definido por ciertas
variables de una base de datos. Dicho lienzo puede ser completado con la
adición de distintas funciones de representación de los datos, como
geom_point() o geom_line() que añaden una nube
de puntos o una línea, respectivamente. Asimismo, funciones cómo
theme() se encargan de cuidar los detalles estéticos del
gráfico, como el tamaño de la letra de los rótulos, el color de fondo o
la presencia o ausencia de ejes. Para más comodidad, en lugar de añadir
la estética para cada gráfico, podemos utilizar la función
theme_set() para definir una estética global que se aplica
a todos los gráficos.
# estética general de los plots
theme_set(theme_minimal() +
theme(axis.title = element_text(size = 15),
axis.line = element_line(color = "black", linewidth = 0.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(color = "black", size = 12),
strip.text.x = element_text(size = 12),
axis.ticks = element_line(color = "black")))# plot de abundancia por sitio
L %>%
# pasamos los nombres de fila a una columna
rownames_to_column(var = "site") %>%
# transformamos la matrix para poder hacer el plot
pivot_longer(cols = -site,
names_to = "species",
values_to = "abundance") %>%
# añadimos el plot
ggplot(aes(x = abundance,
y = reorder(species, abundance, FUN = max),
color = site)) +
geom_point() + # gráfico de puntos
labs(x = "Abundancia",
y = "Organismo",
color = "Sitio") +
theme(text = element_text(size = 18))Queremos entender cómo se distribuyen los organismos en términos de abundancia. A continuación, vamos a ir construyendo nuestro entendimiento de cómo esto sucede, desde las preguntas y análisis más sencillos.
Lo primero de todo: hemos muestreado diez localidades en cinco espacios naturales diferentes. ¿Está la abundancia de los diferentes organismos determinada por el lugar de muestreo?
El paquete mvabund nos permite ajustar modelos para
explicar la abundancia de los organismos. Para esto, es necesaria una
pequeña transformación de la estructure de la matriz de datos
L, que se realiza con la función homónima
mvabund(), integrada en el paquete. Tras esto, somos
capaces de utilizar otra función de mvabund,
manyglm(), para explicar la abundancia de la matriz en
función del lugar de muestreo, explícito como rownames(L)
o, lo que es lo mismo, los nombres de las filas de la matriz L.
Podemos especificar la distribución de los datos para el modelo
utilizando el argumento family; en este caso, lo que mejor
se ajusta es una distribución poisson, pero podéis probar a cambiar esta
distribución (por ejemplo, a family = "negative_binomial")
y ver qué ocurre con la distribución de los residuos.
# creamos un objeto 'mvabund' a partir de la abundancia
spp <- mvabund(L)
# hacemos un modelo para explicar la abundancia en función del sitio
mod <- manyglm(spp ~ rownames(L), family = "poisson") # distr. poisson
# representamos los residuos del modelo
plot(mod)Luego, simplemente lo podemos resolver con un análisis de la varianza:
# análisis de la varianza para ver si el sitio explica la abundancia
(anov <- anova(mod, nBoot = 999))## Time elapsed: 0 hr 0 min 1 sec
## Analysis of Deviance Table
##
## Model: spp ~ rownames(L)
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 9
## rownames(L) 0 9 156.5 0.201
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.
Como podemos ver, el estadístico nos indica que p = 0.201 > 0.05 y, por tanto, podemos decir que el sitio no tiene un efecto significativo a la hora de explicar las abundancias de los organismos.
En segundo lugar, podemos realizar la aproximación más clásica en modelos de distribución de especies: intentar entender qué variables climáticas explican la abundancia de los diferentes organismos. Ciertos modelos más o menos complejos incluso permiten hacer predicciones sobre la distribución potencial de los organismos a través de modelización ambiental.
En nuestro caso, podemos utilizar una función integrada en el paquete
mvabund llamada traitglm(). Esta función
aplica modelos lineales generalizados (GLMs)
para explicar una matriz de datos en función de otra (u otras) matrices
de datos de estructura compatible. Por ejemplo, no podremos explicar una
matriz de 10x20 en función de una matriz de 5x15, pero sí una de 10x20
en función de una de 10x15.
# explicamos la matriz de abundancia en función del ambiente
sdm_env <- traitglm(L, R) # modelo lineal generalizado (glm)## No traits matrix entered, so will fit SDMs with different env response for each spp
# transformar y dataviz
sdm_env$fourth.corner %>%
as_tibble(rownames = "species") %>%
pivot_longer(cols = -species, names_to = "env_var", values_to = "value") %>%
mutate(species = as.factor(str_extract(species, "sp\\d+")),
value = round(value, 5),
env_var = forcats::fct_inorder(env_var),
species = forcats::fct_reorder(species, as.numeric(str_extract(species, "\\d+")))) %>%
# mapa de calor
ggplot(aes(x = env_var, y = fct_rev(species), fill = value)) +
geom_raster() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0) +
labs(x = "Variables ambientales",
y = "Especies",
fill = "Coeficiente") +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))A pesar de que esta aproximación nos aporta gran información sobre la relación de las especies con su ambiente, no deja de ser una forma de analizar qué ocurre con las especies de forma individual. Esto no nos permite alcanzar el nivel de abstracción que deseamos para resolver a nuestra pregunta; para ello, necesitamos incluir los datos de los rasgos funcionales de las especies.
Esto lo podemos hacer con la misma función, traitglm(),
pero añadiendo la matriz Q. En este caso, necesitamos
introducir la transpuesta (t()) de Q porque sino
el modelo no puede computar matrices de diferente tamaño y nos da el
siguiente error:
number of rows of matrices must match (see arg 2).
# añadimos la matriz de rasgos funcionales
f.corner <- traitglm(L, R, t(Q)) # glm de nuevo
# observemos la matriz fourth-corner
f.corner$fourth.corner## env1 env2 env3 env4 env5
## t1 0.005896477 -0.030156079 0.032728336 -0.022556718 -0.026917732
## t2 -0.015514887 0.013106324 -0.012049221 -0.039692052 -0.007347624
## t3 0.020185284 0.024309594 0.049111168 0.001369728 -0.046107386
## t4 -0.004948459 0.012092828 -0.030603234 -0.014503910 0.024181498
## t5 -0.020400040 0.005302025 0.020634272 -0.017760382 -0.009821665
## t6 -0.042718373 0.017132002 -0.056464978 -0.041663403 0.040201806
## t7 0.016447268 -0.032752837 0.066954612 0.021185391 -0.057901879
## t8 0.006364846 -0.025288544 0.051488167 0.018310605 -0.029164103
## t9 0.014388225 -0.047493117 0.024863594 0.010236594 0.014282519
## t10 -0.013240487 0.036835328 -0.015586492 -0.008098365 -0.008877121
## t11 -0.014959676 -0.013503966 0.008225036 0.007433749 -0.012406317
## t12 0.003445801 -0.036153152 0.051730749 0.009830930 -0.019343220
## t13 -0.019283181 -0.024839259 0.048699685 -0.016277801 -0.036393665
## t14 -0.024046745 0.050593395 -0.028186839 -0.009237108 0.013941703
## t15 -0.004714214 -0.001967574 -0.007937043 -0.026548242 0.027757722
## env6 env7 env8 env9 env10
## t1 -0.019726434 -0.029070770 0.015378237 -0.033377981 -0.007832847
## t2 0.001187777 -0.021322501 0.035382373 -0.026892351 0.015762487
## t3 0.028208979 -0.008266815 -0.008361714 0.011775754 0.032692182
## t4 0.004730038 0.018637453 0.007642035 -0.002224654 -0.015310601
## t5 -0.015122491 0.024281807 0.001866518 -0.009335415 -0.019099460
## t6 -0.020126082 0.072155957 0.013309618 0.008538976 -0.026845870
## t7 -0.025845623 -0.002285457 -0.034460053 0.008901942 -0.005117371
## t8 -0.010154626 -0.037706166 -0.010054318 -0.012561046 0.007880418
## t9 0.010102239 -0.080896695 0.030564598 -0.030114349 0.037166281
## t10 0.002662576 0.038930905 -0.002267904 0.016424265 -0.003331159
## t11 -0.033098305 0.024245015 -0.013565228 0.009772219 -0.013308426
## t12 -0.012807564 -0.041765937 0.031867364 -0.015526612 0.037339151
## t13 -0.040957657 0.010185582 -0.017816583 -0.018358099 -0.031138063
## t14 0.001621939 0.028463857 0.031573752 0.001838934 0.005727427
## t15 0.010479725 0.005533823 0.012901708 -0.014480981 -0.010271220
## env11 env12 env13 env14 env15
## t1 0.028755140 0.027739712 -0.001766501 -0.008262744 -0.011906722
## t2 0.018141164 0.050584237 0.034553432 0.011548849 -0.013527046
## t3 0.001882132 0.004668550 0.025864721 -0.003212201 0.030810105
## t4 -0.006084416 0.009707665 -0.011986105 -0.005002110 -0.042549437
## t5 0.004529213 0.019836164 0.011786953 -0.006707977 -0.024753065
## t6 0.006919360 0.036652412 0.013557160 0.022501647 -0.060949812
## t7 0.015657700 -0.018618540 -0.001088288 0.003257411 0.038172591
## t8 0.011025947 -0.015608570 0.010151306 -0.001234242 0.053281915
## t9 0.014327452 -0.009588088 -0.006568942 -0.020472294 0.065077764
## t10 -0.015113117 0.015529544 0.012788965 0.003663808 -0.032685636
## t11 0.013128807 -0.003697762 0.009224238 0.023593957 0.012053497
## t12 0.026922428 0.007788756 0.008441355 -0.009807739 0.055847632
## t13 0.026987546 0.014832604 0.021859639 0.013164277 0.003823548
## t14 0.001589585 0.027440968 0.024061261 0.018156810 -0.030603694
## t15 0.006159755 0.014595450 -0.005079241 -0.012195386 -0.026284344
Y, de esta forma, somos capaces de calcular las interacciones entre los rasgos funcionales de los organismos y el ambiente en el que ocurren. Vamos a representar estos coeficientes gráficamente, no sin antes hacer una rápida comprobación de los residuos del modelo:
# finalmente, podemos representar el mapa de calor
f.corner$fourth.corner %>%
as_tibble(rownames = "traits") %>%
pivot_longer(cols = -traits, names_to = "env_var", values_to = "value") %>%
mutate(value = round(value, 5),
env_var = forcats::fct_inorder(env_var),
traits = forcats::fct_reorder(traits, as.numeric(str_extract(traits, "\\d+")))) %>%
# mapa de calor
ggplot(aes(x = env_var, y = fct_rev(traits), fill = value)) +
geom_raster() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0) +
labs(x = "Variables ambientales",
y = "Rasgos funcionales",
fill = "Coeficiente") +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))¿Qué observamos en esta matriz? Se ve que todas las posibles interacciones han sido otorgadas un determinado valor. ¿Es esto realista? ¿Son todas las interacciones entre rasgos funcionales y variables ambientales importantes?
Lasso (least square shrinkage and selection operator) (Tibshirani 1996) es una forma de regresión que añade una penalización a los coeficientes de mínimos cuadrados, que es el cálculo que se utiliza para encontrar modelos de regresión por comparación de los valores predichos con los observados. Esta penalización puede modificar los coeficientes de las variables del modelo y está controlada por un parámetro, \(\lambda\), que puede tomar valores desde 1 hasta infinito positivo. A medida que \(\lambda\) aumenta, la pendiente de la variable en el modelo mengua. La grandísima ventaja de esto es que lasso es capaz de llevar a cero la pendiente de aquellas variables del modelo que no sean relevantes para explicar los datos, dejando solo aquellas que sí lo sean. Esto nos permite simplificar los modelos y, en general, nuestro trabajo y toma de decisiones. Como recurso, os dejo este vídeo sobre regularización lasso que, aunque en un tono no demasiado serio, consigue explicar bien qué utilidad tiene y cómo funciona. Es ese vídeo oiréis hablar también de la ridge regression, que es muy similar a la regresión lasso pero no es capaz de llevar los coeficientes de las variables poco explicativas a cero.
Lasso está disponible para la aproximación fourth-corner a
través del algoritmo glm1path de la función
traitglm(), como se puede ver a continuación. Sin embargo,
lasso no es un método exclusivo de fourth-corner, sino que se
puede aplicar a cualquier regresión que necesite selección de variables,
como por ejemplo en el análisis de big
data.
# añadimos lasso a través del argumento method = "glm1path"
lasso <- traitglm(L, R, t(Q), method = "glm1path")
# observemos la matriz fourth-corner con lasso
lasso$fourth.corner## env1 env2 env3 env4 env5 env6 env7 env8 env9 env10
## t1 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t2 0 0.000000000 -0.0126001270 0 0 0 0.000000000 0 0 0
## t3 0 0.000000000 0.0007411432 0 0 0 0.000000000 0 0 0
## t4 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t5 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t6 0 0.000000000 0.0000000000 0 0 0 0.009344633 0 0 0
## t7 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t8 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t9 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t10 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t11 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t12 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t13 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## t14 0 0.006611215 0.0000000000 0 0 0 0.000000000 0 0 0
## t15 0 0.000000000 0.0000000000 0 0 0 0.000000000 0 0 0
## env11 env12 env13 env14 env15
## t1 0 0 0 0 0.000000000
## t2 0 0 0 0 0.000000000
## t3 0 0 0 0 0.000000000
## t4 0 0 0 0 0.000000000
## t5 0 0 0 0 0.000000000
## t6 0 0 0 0 0.005031213
## t7 0 0 0 0 0.000000000
## t8 0 0 0 0 0.000000000
## t9 0 0 0 0 0.000000000
## t10 0 0 0 0 0.000000000
## t11 0 0 0 0 0.000000000
## t12 0 0 0 0 0.000000000
## t13 0 0 0 0 0.000000000
## t14 0 0 0 0 0.000000000
## t15 0 0 0 0 0.000000000
# y volvemos a representar el mapa de calor
lasso$fourth.corner %>%
as_tibble(rownames = "traits") %>%
pivot_longer(cols = -traits, names_to = "env_var", values_to = "value") %>%
mutate(value = round(value, 5),
env_var = forcats::fct_inorder(env_var),
traits = forcats::fct_reorder(traits, as.numeric(str_extract(traits, "\\d+")))) %>%
# mapa de calor
ggplot(aes(x = env_var, y = fct_rev(traits), fill = value)) +
geom_raster() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0) +
labs(x = "Variables ambientales",
y = "Rasgos funcionales",
fill = "Coeficiente") +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))La semilla que sembramos al principio, 5558, está seleccionada para generar unos datos que producen estos resultados concretos. ¿Qué ocurre si cambiamos la semilla? Probad a rehacer los análisis con semillas que tomen los siguientes valores: 7120, 9937, 1840, 5647 y 216.
En esta segunda parte vamos a dejar de lado los datos dummy
y, en su lugar, vamos a analizar una base de datos real
de macroinvertebrados de agua dulce. Se trata de 76 familias de
organismos, desde efímeras hasta gammáridos. La lista completa se puede
encontrar usando el siguiente comando: colnames(L), siempre
después de haber cargado la matriz L correspondiente (siguiente paso).
Asimismo, los datos proceden de cinco diferentes espacios naturales de
los Montes Vascos, a saber: AIZK =
Aizkorri (Gipuzkoa), ARAL =
Aralar (Gipuzkoa), ARTI =
Artikutza (Nafarroa), GORB =
Gorbeia (Bizkaia/Araba), e IZKI =
Izki (Araba).
Cargamos la matriz L con el siguiente comando:
Cargamos ahora la matriz Q de rasgos funcionales por familia:
En la tabla siguiente se muestra la descripción de los rasgos
funcionales de los organismos. A tener en cuenta: los valores de los
rasgos se corresponden con los valores mostrados como columnas. Esto
quiere decir que, por ejemplo, si la variable 1 (var1)
tiene un valor de 1.33, ese organismo puede alcanzar un tamaño máximo de
entre 0.25cm y 5cm, pero es más probable que su tamaño sea de 0.25cm o
menor (porque 1.33 está más cerca de 1 que de 2). A esto se le llama fuzzy-coding
los datos de rasgos funcionales.
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|
| Tamaño potencial máximo del organismo | |||||||
| var1 | X0.25cm | X0.250.5cm | X0.51cm | X12cm | X24cm | X48cm | X8cm |
| Distribución longitudinal | |||||||
| var13 | crenon | epirithron | metarithron | hyporithron | epipotamon | metapotamon | estuary |
| Altitud en la que se encuentra | |||||||
| var14 | lowlands | piedmont.level | alpine.level | ||||
| Velocidad de corriente (preferencia) | |||||||
| var17 | null | slow.less.25.cms | medium.25.50.cmlesss | fast.more.50.cms | |||
| Status trófico (preferencia) | |||||||
| var18 | oligotrophic | mesotrophic | eutrophic | ||||
| Saprobicidad (preferencia) | |||||||
| var21 | xenosaprobic | oligosaprobic | b.mesosaprobic | a.mesosaprobic | polysaprobic | ||
| pH del medio (preferencia) | |||||||
| var22 | less.4 | more.4.4.5 | more.4.5.5 | more.5.5.5 | more.5.5.6 | more.6 | |
Y, finalmente, cargamos la matriz R con las variables ambientales:
En la siguiente tabla se proporciona información sobre el significado de las variables ambientales medidas en los espacios naturales muestreados.
| Variable | Descripción |
|---|---|
| pH | Valores de pH en el medio |
| Conductivity_.µS.cm. | Conductividad (microS/cm) |
| O2_sat | Saturación de oxígeno (%) |
| O2_conc | Concentración de oxígeno (mg/L) |
| Cl_.mg.L._avg | Concentración de cloruro (mg/L) |
| SO4_.mg.L._avg | Concentración de sulfato (mg/L) |
| NO3_.mg.L._avg | Concentración de nitrato (mg/L) |
| N.NH4_.mg.L._avg | Concentración de amonio (mg/L) |
| DIN.N.mg.L. | Concentración de nitrógeno (mg/L) |
| P.PO4_.µg.L. | Concentración de fosfato (mg/L) |
| tree_cover_. | Cobertura arbórea (%) |
| BARE_SOIL_._19_23 | Cobertura de suelo descubierto (%) |
| GRASSLAND_._19_23 | Cobertura de praderas (%) |
| SHRUBLAND_._19_23 | Cobertura arbustiva (%) |
| NAT_FOREST_._19_23 | Cobertura de bosque natural (%) |
| PLANT_FOREST_._19_23 | Cobertura de plantaciones (%) |
| Catchment_area_km2 | Área de muestreo (km2) |
El objetivo es calcular las interacciones significativas entre los
rasgos funcionales y las variables ambientales que nos ayuden a explicar
la distribución de la abundancia de los diferentes organismos. Para
ello, comenzaremos visualizando los datos y, después, utilizaremos la
función traitglm() del paquete mvabund tal y
como vimos anteriormente: primero, sin penalización lasso y, segundo,
con ella para seleccionar interacciones.
Vamos a hacer un par de gráficos que nos ayuden a entender mejor qué datos tenemos.
# plot de abundancia por sitio
L %>%
# rownames (sites) to column to plot
rownames_to_column(var = "site") %>%
# transformamos la matrix para poder hacer el plot
pivot_longer(cols = -site,
names_to = "species",
values_to = "abundance") %>%
# añadimos el plot
ggplot(aes(x = abundance,
y = reorder(species, abundance, FUN = max),
color = site)) +
geom_point() + # gráfico de puntos
labs(x = "Abundancia",
y = "Familias",
color = "Sitio") +
theme(text = element_text(size = 18))¿Se parece esta distribución de las abundancias a la que teníamos anteriormente para los datos dummy?
Vamos ahora directamente a resolver la pregunta que nos interesa,
para lo cual utilizaremos las tres matrices como argumentos de
traitglm(), y representamos gráficamente la matriz de
interacciones seleccionada. Lo hacemos para ambos modelos, sin lasso y
con lasso (recordamos: method = glm1path) y los
representamos.
Primero sin lasso:
f.corner <- traitglm(L, R, Q)
f.corner$fourth.corner %>%
as_tibble(rownames = "traits") %>%
pivot_longer(cols = -traits, names_to = "env_var", values_to = "value") %>%
mutate(value = round(value, 5),
env_var = factor(env_var, levels = sort(unique(env_var))),
traits = factor(traits, levels = sort(unique(traits), decreasing = TRUE))) %>%
# mapa de calor
ggplot(aes(x = env_var, y = traits, fill = value)) +
geom_raster() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0) +
labs(title = "Modelo independiente del espacio natural",
x = "Var. ambientales",
y = "Rasgos func.",
fill = "Coef.") +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))Y después con lasso:
lasso <- traitglm(L, R, Q, method = "glm1path")
lasso$fourth.corner %>%
as_tibble(rownames = "traits") %>%
pivot_longer(cols = -traits, names_to = "env_var", values_to = "value") %>%
mutate(value = round(value, 5),
env_var = factor(env_var, levels = sort(unique(env_var))),
traits = factor(traits, levels = sort(unique(traits), decreasing = TRUE))) %>%
# mapa de calor
ggplot(aes(x = env_var, y = traits, fill = value)) +
geom_raster() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0) +
labs(title = "Modelo independiente del espacio natural",
x = "Var. ambientales",
y = "Rasgos func.",
fill = "Coef.") +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))¿Qué está pasando aquí? ¿Por qué parece que hay menos interacciones seleccionadas para el gráfico sin lasso (izquierda) que para el gráfico con lasso (derecha)? Vamos a ver las matrices de interacciones:
## pH Conductivity_.µS.cm. O2_sat O2_conc Cl_.mg.L._avg
## var1 0.01315904 0.08485679 -0.04385779 -0.09429455 0.159752592
## var13 0.03724744 0.07064226 -0.12026136 0.10177320 -0.130719892
## var14 0.01182520 -0.20460362 -0.06799693 0.05762894 -0.068914605
## var17 -0.22817063 0.03082527 0.36298936 -0.05095678 -0.308052877
## var18 -0.04010953 0.09273085 0.12981674 -0.03694424 -0.066904964
## var21 -0.15739570 -0.01029542 0.05778800 -0.14558510 0.037993393
## var22 -0.14464822 0.08830754 0.09045460 -0.02330884 -0.003115431
## SO4_.mg.L._avg NO3_.mg.L._avg N.NH4_.mg.L._avg DIN.N.mg.L. P.PO4_.µg.L.
## var1 -0.13424086 -0.07621326 0.038806233 -0.07479167 -0.12263507
## var13 -0.01803699 0.13607971 0.210911512 0.14646150 0.01439262
## var14 0.10573982 -0.03919476 0.026474914 -0.03816324 0.02472374
## var17 -0.05487787 0.08396887 0.020640364 0.08532485 -0.14851299
## var18 -0.12750873 -0.08065156 0.001747011 -0.08096013 -0.11096515
## var21 0.05871555 -0.03579921 -0.045194207 -0.03805644 0.04864642
## var22 -0.08757983 0.07776019 -0.260557410 0.06612103 0.10116822
## tree_cover_. BARE_SOIL_._19_23 GRASSLAND_._19_23 SHRUBLAND_._19_23
## var1 0.079114951 -0.11543830 -0.1971582 -0.2500568
## var13 0.012928469 0.29542957 0.9264349 1.7104471
## var14 -0.005695622 -0.06749943 -0.3468541 -0.8438421
## var17 0.041052274 -0.01778968 0.2630383 0.9119774
## var18 0.064415597 0.02771697 0.3715773 0.7648960
## var21 -0.060071252 -0.18598408 -0.5006729 -1.3722988
## var22 -0.122324576 -0.06634006 -0.1087394 -0.4405074
## NAT_FOREST_._19_23 PLANT_FOREST_._19_23 Catchment_area_km2
## var1 -0.4186024 -0.1426231 -0.01484885
## var13 2.4691483 1.7927658 0.06448584
## var14 -1.1708913 -0.8676686 -0.15495297
## var17 1.4136504 0.9392763 0.22122586
## var18 1.3579632 1.0860101 0.05005877
## var21 -1.8517342 -1.3649009 0.06347529
## var22 -0.7797123 -0.4636749 0.01265808
## pH Conductivity_.µS.cm. O2_sat O2_conc Cl_.mg.L._avg
## var1 0.00000000 0.04532127 -0.04740023 -0.040929013 0.13017623
## var13 0.00000000 0.00000000 -0.03418713 0.000000000 -0.02853290
## var14 0.00000000 -0.07377945 -0.02252036 0.000000000 -0.07853568
## var17 -0.06788648 -0.13447846 0.15216154 0.015706817 -0.04396393
## var18 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000
## var21 -0.01216116 0.02524483 0.00000000 -0.054184873 0.02458062
## var22 -0.02059854 0.00000000 0.06917813 0.003542237 0.00000000
## SO4_.mg.L._avg NO3_.mg.L._avg N.NH4_.mg.L._avg DIN.N.mg.L. P.PO4_.µg.L.
## var1 0 -0.09755478 0.016320134 0.00000000 -0.069472661
## var13 0 0.00000000 0.076288129 0.07944128 0.020166768
## var14 0 0.00000000 0.001604626 0.00000000 0.000000000
## var17 0 0.01974271 0.000000000 0.00000000 -0.004395314
## var18 0 0.00000000 0.000000000 -0.10166980 0.000000000
## var21 0 0.00000000 0.000000000 0.00000000 0.000000000
## var22 0 0.06100831 -0.153173000 0.00000000 0.022676177
## tree_cover_. BARE_SOIL_._19_23 GRASSLAND_._19_23 SHRUBLAND_._19_23
## var1 0.03139716 -0.03484952 -0.08209893 0.000000000
## var13 0.00000000 0.00000000 0.16596512 0.000000000
## var14 0.01500141 0.00000000 0.00000000 0.000000000
## var17 0.00000000 -0.05502395 -0.03508282 0.002630367
## var18 0.00000000 -0.01780711 0.00000000 -0.056330545
## var21 -0.02603755 0.00000000 0.00000000 -0.069400967
## var22 -0.03645761 0.00000000 0.03377143 0.000000000
## NAT_FOREST_._19_23 PLANT_FOREST_._19_23 Catchment_area_km2
## var1 -0.021488948 0.02589767 -0.009847158
## var13 0.000000000 0.00000000 0.093095220
## var14 0.000000000 0.00000000 -0.125846106
## var17 0.000000000 -0.01222346 0.136464856
## var18 0.001169351 0.00000000 0.000000000
## var21 0.000000000 0.00000000 0.000000000
## var22 -0.037800510 0.00000000 0.034095625
Como vemos, hemos de ir con ojo con las representaciones visuales, ya que una escala de colores puede darnos una idea errónea de los resultados.
El objetivo de este ejercicio final es resolver la siguiente pregunta: ¿qué ocurre con las interacciones rasgo-ambiente si analizamos cada espacio natural por separado?
Pista: primero tenemos que seleccionar solamente las filas que correspondan al espacio natural que nos interese. Esto lo podemos hacer con la función
grepl()del siguiente modo:X[grepl("CODE", rownames(X)), ], siendoXla matriz que corresponda y siendoCODEel código del espacio natural deseado. El resto del proceso es idéntico a lo explicado anteriormente:traitglm()con las tres matrices y el argumento para usar lasso.
Aquí los resultados:
Basque Center for Climate Change (BC3); rodrigo.granjel@bc3research.org↩︎